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We report an extensive numerical study of a charged colloidal system with competing short-range 
depletion attraction and long-range electrostatic repulsion. By analizing the cluster properties, we 
identify two distinct regions in the phase diagram: a state composed of stable finite-size clusters, 
whose relative interactions are dominated by long-range repulsion, and a percolating network. Both 
states are found to dynamically arrest at low temperatures, providing evidence of the existence of 
two distinct non-ergodic states in these systems: a Wigner glass of clusters and a gel. 

I. INTRODUCTION 

The origin of low-density non-ergodic states in colloidal systems is a matter of continuous debate and ongoing 
research [H [21 [3l HJ [5]. Several different mechanisms may concur in the formation of these arrested states, depending 
on the relative ratio between the thermal and binding energy and on the shape and symmetry of the interaction 
potential. When particles interact via excluded volume interactions complemented by a spherically symmetric, at- 
tractive potential, it has been shown that the formation of a gel structure takes place concurrently with a spinodal 
decomposition process. Gelation results from an arrested phase separation [6. . A different scenario occurs when col- 
loidal particles have a residual electrostatic charge which builds up an additional long-range repulsion in the effective 
colloid-colloid interaction. This term is often modeled as a Yukawa potential with Debye screening length £ to take 
into account the presence of the solvent and counterions[7]. In apolar solutions or under low salt conditions, when 
particles are sub-micron sized, £ can become comparable to the particle dimension. This long-range repulsive term 
can coexist with a short-range attraction (which can be induced for example via depletion interactions), generat- 
ing a competition between aggregation driven by the attractive part of the potential and the stabilizing role of the 
repulsion, which may ultimately suppress the macroscopic phase separation. Indeed, it has been shown that the 
addition of a long-range repulsion of moderate strength can shift to larger attraction strengths (or lower tempera- 
tures) the phase separationjHJ [S], eventually inhibiting it [TDJ [HI [HI E3]- In this case a microphase separation into 
clusters [HI [T51 [TU [TTJ UH] of a preferred cluster size and shape takes place, depending on repulsion parameters [TT)]. 
When repulsion is moderately short-ranged, i.e. £/er < 0.5 with a being the diameter of the colloidal particle, it was 
observed both in experiments and in simulations that elongated clusters are formed at low enough temperature T 
[2TJ1 [21] • The repulsion between such clusters is relatively weak so that they tend to form at low T quasi-ordered 
columnar structures [22 23J. At large enough packing fraction, the clusters are found to merge into a percolating 
network (20] ETJ . This network of clusters exists at low enough T and undergoes dynamical arrest, so that a gel state 
can be properly identified [I]. 

For cases where the repulsion term is considerably longer-ranged, i.e. £/<r > 1, simulations at low enough colloidal 
densities [Hj have reported the presence of a Wigner glass of clusters. This corresponds to a disordered state of 
polydisperse clusters (due to finite T) which do not percolate and are actually arrested due to the long-range repulsion, 
in analogy with the Wigner glass reported by Chaikin and coworkers [251 Erj] for charged colloidal particles under very 
dilute conditions, stabilized by the Coulomb repulsion. Recently, a comparison between theory and simulations [2"T] 
of Yukawa particles has shown that the ideal Mode Coupling Theory (MCT) provides a quite accurate description 
of the formation of a particle Wigner glass. The MCT predictions for Yukawa particles have also been exploited for 
interpreting arrest into a Wigner glasses of clusters, in systems with competing interactions [24, 28J. Indeed, once 
clusters are assumed spherical and monodisperse (in size), the effective cluster-cluster interactions can be modeled in 
terms of a Yukawa potential, with the same screening length as the one acting between single particles but with a 
renormalized amplitude 24J . 

These earlier works call for additional investigations, in order to further question the existence and the stability 
of a Wigner glass of clusters, as well as a deeper understanding of cluster-cluster interactions. To this end, it is also 
relevant to mention a recent simulation study |29j where clusters were observed to arrest, at not-too-low density, by 
percolation rather than by repulsion. However, differences in the simulation protocol of this work are present with 
respect to that used in [23], in particular history of quench and quench rate, as well as a shorter cut-off distance for 
the long-range repulsion potential. Hence, a more comprehensive study of these models in a wide region of packing 
fraction (f> and temperature T, fully accounting for the long range nature of the repulsive interactions, is needed. 
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In this work we report an extensive simulation study which aims at elucidating in detail the phase diagram of 
colloidal systems interacting with both short-range attraction and long-range repulsion, encompassing states well 
below, at the crossing and well above the percolation line at low enough T. We characterize in detail the properties 
of the aggregates that are formed, both when they exist in finite size objects (clusters) and when they merge into a 
percolating network. We analyze both particle-particle and cluster-cluster correlations to show that, at low enough 
packing fractions, the system self-organizes into stable, long-living clusters which do not percolate and do not form 
ordered structures. The interactions between these clusters can be characterized, for low and intermediate <f>, in terms 
of renormalized long-range repulsive interactions of Yukawa form, as previously hypothesized [24 . Most importantly, 
despite the high polydispersity and non-sphericity of the clusters, we find that the screening length of cluster-cluster 
interactions remains unchanged with respect to that of particle-particle interactions, while the repulsion amplitude is 
found to increase with particle density. At larger tfi, when the shape of the clusters starts to significantly deviate from 
the spherical one and increased packing leads to the arising of branching events, a crossover takes place, ultimately 
leading to the percolation of the clusters associated to a gel transition. We also monitor the dynamics of the two 
regions, calling for the existence of two distinct non-ergodic states in this low-density, low-T part of the phase diagram. 



II. DESCRIPTION OF THE SIMULATION PROTOCOL AND ANALYSIS 



We consider a system of N = 1000 particles of unit mass m interacting with a total pair potential composed of 
a short-ranged attractive part, modeled for convenience as a generalized Lennard- Jones potential [30j with exponent 
a = 100, and of a long-range screened electrostatic repulsion, modeled as a Yukawa term, 

"M -*[(?)""-(?) 

where A is the amplitude of repulsion and e, the depth of the attractive part, is chosen as the energy unit. Time 
is measured in units of yjma 2 Je. The choice of a — 100, also studied in [Tj3 [H], ensures a very short-range 
attraction, corresponding to bond formation only within the first neighbour shell. According to the extended law 
of corresponding states for spherical short-ranged attractions [3T], such choice is generic for any width and shape of 
the attractive potential[6]. Note that also the cluster ground-state properties have been shown to be invariant for 
a > 18 [19]. We fix the parameters of the Yukawa potential to £ = 2er and A — 0.2e. For such values a microphase 
separation into clusters exists [HO HJJ. To correctly take into account the long-range nature of the interactions we 
solve the equations of motion using Ewald summation |32J. Indeed, we have compared results obtained in this way 
with those based on the use of a finite (although large) cut-off. For our choice of parameters, an underestimate of 
about 10% for the repulsive potential energy is provided by the use of a finite cutoff at 8£, in agreement with previous 
studies |33J . Hence, despite the significantly increased computational cost, it appears crucial for the case under study 
to treat the long-range repulsive term with Ewald sum, in order to discriminate cases where clusters truly form a 
disconnected or a percolating state and to address satisfactorily the nature of the arrested state (s) for the chosen 
value of £. 

To model the motion of colloidal particles in a solvent, we use Brownian dynamics simulations with time step 
St = 0.005 and bare diffusion coefficient Dq = 0.005. For this choice, the crossover from ballistic to diffusive regime, 
for isolated particles, takes place for t ~ 10. We do not treat explicitly the effect of the solvent, i.e. hydrodynamic 
interactions are neglected, but we do not expect these to provide significant changes to the long-time structures that 
we observe [S3]. We neglect any change in the electrostatic parameters, in particular of £, with increasing colloid 
packing fraction <p. For particles of the order of /ito, this variation is expected not to be significant with respect to 
the particle diameter |35j. 

The system was initially prepared at several densities and high T, and later it was slowly equilibrated to successively 
lower T. Our aim is to study the approach (from the equilibrium side) to a dynamic arrest transition, rather than 
a rapid quench inside the region where arrest is observed. We study 9 isochores and several T in the low-0 region 
up to = 0.20. Equilibration was carried out in a NVT ensemble, followed by production runs for data collection 
and analysis. For very low T (depending on the studied 0), a true equilibration is no longer possible. Monitoring for 
example the energy per particle, the system at first shows a robust decrease towards an apparent equilibrium state, 
but then the system starts to display a very slow (logarithmic in time) energy drift, typical of an approach to dynamic 
arrest and trapping in a metastable state. 

The connectivity properties of the system have been monitored by studying the intcrparticlc bonding. Two particles 
are considered bonded when the distance between them is smaller than the position of the local maximum in the 
interaction potential, i.e. r < — 1.072c. When a state point has reached equilibrium, we collect several independent 
realizations of the system and calculate the distribution n(s) of clusters of size s over time. Moreover, we check whether 
the largest cluster that we find in the configuration has spanned the whole box at least in one direction. If at least 
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50% of the considered realizations contains a spanning cluster, the state point is a classified as percolating. For state 
points which do not percolate at the studied T, we monitor the cluster properties in time, finding that the clusters 
have a finite, but large lifetime. If the cluster size distribution shows a maximum for a finite value of s, indicating 
that a preferred size for the clusters exists, we classify this state point in the cluster phase region. 



III. FINITE SIZE CLUSTERS AND PERCOLATION: DISORDERED, METASTABLE STATES 

In Fig. [l]we report the phase diagram of the system in the studied (4>,T) region, highlighting the loci of points 
where (i) a cluster phase is observed and (ii) the system percolates. We also report a (dashed) line delimiting the 
region where a true equilibration cannot be reached, meaning that, for all investigated state points below the line, the 
potential energy keeps displaying a very slow, logarithmic drift for the entire duration of the simulation. This region 
largely belongs to the portion of the phase diagram where stable clusters and percolation are present. 




FIG. 1: Phase diagram of the studied system. Open symbols indicate fluid states (no cluster phase and non-percolating). Filled 
diamonds indicate points within the cluster phase region (i.e. showing a maximum at finite s for n(s)) and filled squares refer 
to percolating state-points. Full lines, defining the cluster region and the percolation region, and the dashed line, defining the 
region (i.e. the state points below the line) where the energy starts to display a slow logarithmic drift, are drawn as guides to 
the eye. At low enough <j> (i.e. (f> < 0.02) the cluster size distribution is peaked at s = 1. Inset: potential energy per particle 
e p versus packing fraction. A minimum is found at very low <f>; e v becomes positive for cf> > 0.125, i.e. close to the location of 
the percolation transition. 



We notice that in the absence of repulsion, the critical temperature for the short-ranged attractive potential studied 
here is found at T c ~ 0.24[2I], while the critical packing fraction is found to be approximately <fi c ~ 0.27, following the 
behavior for spherical attractive potentials in the short-range limit |36j. Hence, the addition of a long-range repulsion 
suppresses the tendency to phase separate macroscopically, at least in the T-window that we have studied, i.e. for 
T > 0.05. It becomes thus possible to reach lower temperatures in one-phase condition. At these low T, the lifetime 
of the inter-particle bond increases significantly and clusters behave as effective long-lived aggregates. 

Fig. [I] shows that for cf> < 0.125 the system never percolates at all studied T. For very small cf> (e.g. 4> = 0.005) 
the system organizes into very small clusters, and the maximum of n(s) is located to s = 1 at all studied T. We 
can consider this very low-</> region to be essentially in a 'monomeric' state. As <f> is increased, a maximum in the 
cluster distribution n(s) develops for T < 0.15, and we can identify a stable cluster phase, according to our definition 
discussed above. In this region, clusters can be monitored in time and we observe that they continuously exchange 
particles with each other, breaking and reforming bonds, while on average existing at all times in different realizations. 
For T < 0.1, clusters become more and more long-lived, since the time necessary to break a bond becomes comparable 
to the simulation time. Similarly, the encounter between two different clusters is quite rare, due to the presence of 
the long-range repulsion. Essentially, clusters become frozen. Finally, for <f> > 0.14 we observe percolating states. 
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Hence, at low T or equivalently for high attraction strengths, the system displays a transition from monomers to a 
stable finite-cluster phase and eventually to a particle network. These findings are in very good agreement with recent 
confocal microscopy experiments carried out for a charged colloidal suspension with depletion interactions under low 
screening conditions [37], consistent with the case discussed in the present work. 

In the inset of Fig. [l] we report the dependence of the potential energy per particle e p , on </> for T = 0.1. Data 
are collected at the end of the long simulation runs (due to the presence of the slow logarithmic aging). We notice 
that the energy has a minimum for 0.02 < (f> < 0.04 and then becomes positive for (f> > 0.125, i.e. in the region 
where percolation is observed. This points to the fact that only at very low <f> particle clustering is dominated by the 
attractive interactions, while the interplay with repulsion becomes important already for 4> ^ 0.04. The increase in 
the energy arises from the large number of neighbours within the repulsive range of the interaction, which overcomes 
the gain associated to bonding. 




FIG. 2: Snaphots of the simulated system for serveral different 4> at T = 0.10. Different clusters are drawn in different colors 
according to particle number in the cluster. For 4> < 0.125 (top row), the system forms at all times several distinct clusters 
which are not connected and do not span the simulation box; while for <fi > 0.14 (bottom row) a percolating cluster exists at 
all times. 

To visualize the cluster formation and to highlight the fact that the system at low </> neither percolates, nor forms 
ordered structures, we report in Fig. [2] some snaphots of the system for different <fi along the isotherm T = 0.10. 
Different clusters are drawn in different colors, according to particle number. Although exchange of particles between 
clusters as well as cluster branching and breaking processes are sometimes observed, a similar picture of clusters to 
that reported in Fig. [2]exists at all times. The emergence of a percolating cluster arises, for the reported temperature, 
at <j> = 0.14, and rapidly involves the majority of particles as <fi grows further. We notice that in the cluster phase 
region, clusters are polydisperse, both in size and in shape. While at low <p they are always rather spherical, they tend 
to become more and more elongated with increasing packing fraction. This phenomenon was observed in previous 
simulation 21J and experimental [20 studies for a system with a much shorter screening length. A study of the ground 
state properties of isolated clusters [19] has shown that, while in the short screening length case quasi one-dimensional 
cluster growth is energetically favoured, giving rise to the peculiar Bernal spirals [2TJI 121], for the present study the 
expected ground state cluster structure is much more spherical, although with some degree of anisotropy. Comparing 
the ground state structure of clusters of a certain size, shown in Fig. 5 of Ref.JTU] for the two cases, it is evident that 
the average number of nearest neighbours is dramatically different: while in the Bernal spiral, particles have always 
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exactly 6 neighbours (since there is no difference between bulk and surface), in our case particles in the interior of the 
cluster have a coordination close to 12 neighbours, while those on the surface a much smaller one (close to 6). Hence, 
the resulting average coordination is close to 8. 

To gain a better understanding of the cluster shape and local structure for the present case, we report in Fig. [3] the 
(average) distribution of nearest neighbours P(n) for all studied <f> for T = 0.1, and the average number of neighbours 
< n > in the associated inset, which is compared to the ground state predictions. While, as expected, the number 
of nearest neighbours at first increases with 0, above percolation it roughly stops evolving, so that < n > does not 
grow much above 6. Indeed, the distribution tends to saturate and to remain always peaked around 6 neighbours, 
despite a gradual increase in the number of particles with large number of neighbours. Monitoring the evolution 
of < n > with T (not shown), we do not observe the presence of non-monotonic effects, which were found in the 
low-screening length case[2T], suggesting that the long-range repulsion never facilitates the formation of compact 
structures. Comparing with the isolated cluster study of Ref ■ [15] . we can argue that the interactions between clusters 
act against the compaction of the clusters and favor the formation of elongated structures. Indeed, if we imagine 
cluster-cluster interactions to be ruled by the long-range repulsion (as we will prove below), clusters will occupy 
much more efficiently the space and, at the same time, reduce the total potential energy due to increased average 
cluster-cluster distance by growing in an elongated rather than spherical manner. 




FIG. 3: Distribution of nearest neighbours P(n) for all studied <j> and T = 0.1. Inset: dependence on <j> of the average number 
of neighbours < n > and comparison with isolated clusters ground state predictions from |19| . 

Next, we examine the behaviour of the radius of gyration of clusters of size s, defined as R g (s) = 

1/2 

jT72"([Si=i( r i — R-cm) 2 ] )j where are the coordinates of particle i, Hcai is the cluster center of mass and 
the average is performed over all particles of size s. We report R g (s) in Fig. [4] for all studied <fi at T — 0.1. Percolating 
clusters are not included in the analysis. 

We notice that, in agreement with data in Ref.|2T], the dependence of R g on cluster size below the percolation 
threshold does not show a dependence on <f>, although larger cj> values allow sampling of larger cluster sizes. However, 
as <f> grows above percolation, clusters gradually become slightly more compact, displaying a smaller R g at comparable 
sizes. We can try to identify a growth law for R g with size. For s < 10, clusters are compact objects, as expected. 
Defining a fractal dimension for clusters df as the power law exponent of R g ~ s 1 ' d l , we find that for s > 10, 
clusters lose their compactness and a quasi-linear growth is observed, roughly compatible with an exponent 1.3. This 
is close to what observed for the Bernal spiral case. The figure also shows the corresponding data for the ground state 
configuration of isolated clusters, taken from Ref. [19 . In comparison with the isolated clusters ground states, the 
clusters observed in the simulations are significantly less compact, suggesting that interactions between clusters and 
the non-negligible role of entropy induce a one-dimensional growth already at smaller sizes. Hence, the system always 
remains very far from isolated clusters ground state predictions. Probably this results from a combination of effects: 
on one hand, the increased bond and cluster lifetime does not allow an effective restructuring of the clusters towards 
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FIG. 4: (a) Radius of gyration of finite-size clusters for all studied <f> at T = 0.1. Also shown are the ground state calculations 
(open circles) for isolated clusters taken from [T§]. Dashed lines indicate power-law behavior s 1 ^ df , namely spherical (df = 3) 
for small sizes and quasi-linear (df = 1.3) for large sizes. 



their preferred configuration, on the other hand entropy and cluster-cluster interactions act against spherical growth. 
In summary, we find the system existing in highly disordered, metastable states, both above and below percolation. 
This is very clear looking at the Movies, attached for several studied cj> at T = 0.1 as supporting material to the present 
work [31]. These results point to the existence of two distinct metastable disordered states, that we are tempted to 
identify as a Wigner glass of clusters below percolation and a gel state above percolation. To verify whether these two 
states are kinetically arrested, we will examine the dynamics of the system in the following. However, we can already 
report that an ordered low-T state is never reached during the course of our simulations at any studied <fi, despite 
the slow equilibration path that we have used, probably due to the large screening length used. Indeed, a recent 
study[23 where a much shorter screening length was used (although the potential studied there was not Yukawa-like 
but exponential) reported the formation of a columnar phase during the course of the simulation: the tendency to 
order was increased in that case by the reduced repulsive barrier, facilitating formation and breaking of bonds. 

Finally, we report the cluster size distribution along the same isotherm, for all studied <fi, in Figure [5j We observe 
a number of relevant phenomena: (i) the emergence of a clear peak between 0.02 < (f> < 0.125, delimiting the cluster 
phase region for this T as we have defined it and whose boundaries have been represented in Fig. [I] The peak arises 
at successively larger sizes as expected; (ii) the gradual disappearance of monomers (s — 1), which from dominant 
at low <f), in what we have called the monomeric phase, become absent. Not only monomers disappear, but also 
small clusters of increasing sizes gradually disappear, and the system exists in the form of finite-size clusters; (iii) for 
increasing 0, numerical noise becomes important, due to the small number of finite clusters, but, within such noise, 
the cluster distribution is consistent with the random percolation power-law prediction, i.e. n(s) ~ s~ 2 - 2 . 

We conclude this paragraph by observing that, in the present study, we do not find any signature of reentrant 
percolation as it was found in the short-range repulsion case |21j . In the latter case, at high enough <j>, percolation 
was observed at first for high T, due to random aggregation of particles, then a restructuring into the preferred shape 
was observed, giving rise to a non-percolating regime for intermediate temperatures, and finally to a new random 
percolation of the spiral- like clusters at low T. Here, in the investigated region (up to <j> = 0.20), we do not observe 
such behavior. 



IV. A MORE CAREFUL LOOK AT THE CLUSTERS: INTRA-CLUSTER AND INTER-CLUSTER 

PROPERTIES 



We have seen so far that, at low </>, the system remains organized into several clusters for all studied temperatures. 
We have tried also deeper quenches (e.g. T — 0.01), where the system remains far-from-equilibrium, and we never 
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FIG. 5: Cluster size distribution for all studied (f> at T — 0.1; only non-percolating clusters are included. The dashed line is 
the prediction for random percolation n(s) ~ s -2 ' 2 . 

observe the coalescence of such clusters due to the long-range repulsion. We have also seen that the clusters are 
polydisperse, even though a preferential size emerges as a peak in the n(s), and that no long-range order is present, 
i.e. we do not observe the presence of any columnar or lamellar phase. However we notice that, at low T, a certain 
degree of order develops inside the clusters. 

By looking at the particle-particle radial distribution function g(r), shown in Fig. [6] for = 0.08, clear sharp peaks 
arise with decreasing T. After a very large maximum at contact (highlighted in the inset), enhanced correlations 
are found in correspondence of specific discrete distances indicating locally preferred geometries (e.g. triangular, 
tetrahedral, linear order, etc.). However, liquid-like disorder is retained after the 2nd peak of the g(r). This behaviour 
is invariant whether we consider states below and above percolation, as also shown in the figure. A parallel analysis 
of the structure factor (not shown) confirms these observations. This indicates that, at low T, particles inside the 
clusters tend to occupy preferential ordered positions, thus providing quasi-crystalline character to the inside of the 
clusters. However, this order is lost already after the second neighbour shell, both in the cluster region and in the 
percolating one. 

After the analysis of particle-particle correlations, it is interesting to consider cluster-cluster correlations to get 
an idea of what mechanisms regulate the interactions between different clusters and why they do not percolate at 
sufficiently low </>. In a previous work [23], it was hyphothesized that cluster interactions may be of renormalized 
Yukawa form, maintaining the same screening length as the underlying particle-particle interactions and with an 
increased repulsion strength with increasing cluster size. This hyphothesis was based on the assumption of spherical 
and monodisperse clusters, an assumption which is not strictly verified, as we have already discussed, within the 
present system. To quantify the inter-cluster interactions, we adopt the following strategy. We consider a low enough 
T where the cluster phase extends over a large <fi region. In each configuration, for each cluster we calculate its 
center of mass coordinate. Then, independently of cluster size and shape, we calculate the pair distribution function 
ffCM 4 ( r ) between centers of mass of different clusters. Results for T — 0.1 and various studied <f> are reported in 
Fig. [7] It is clear that the centers of clusters are found in liquid- like configurations, with increasing correlations as 
(f> increases up to roughly <fi = 0.08. Next we compare the numerical 5 C 7M*( r ) with theoretical predictions obtained 
solving numerically the Ornstein-Zernike (OZ) equation [7] with the hypernetted chain (HNC) closure. We checked 
that results are independent on the chosen closure. As suggested in Ref. [24] , to model the cluster-cluster interaction, 
we select a pure Yukawa potential (no hard-core), with the same screening length (£ = 2.0) of the particle-particle 
repulsive interaction. We leave the amplitude A e ff of the Yukawa potential as the only fitting parameter. Indeed, 
the number density of the clusters is read directly from the simulation data. The best-fit curves are reported in Fig. [7] 
and show a remarkable agreement with the data extracted from the simulations. The behavior of the cluster number 
density and of the renormalized amplitude are reported in the insets. A e ff grows approximately linearly with </>, 
starting from the particle-particle value at <fi —> 0, while the number of clusters progressively shrinks. 
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FIG. 6: Radial distribution function g(r) for particle-particle correlations along the isochore (j> = 0.08 ( within the cluster phase 
region) and several studied T. Also, results for cf> — 0.16 (within the percolating region) at low T (dashed line) are reported to 
show the invariance of the locally preferred positions and liquid-like disorder at large distances in the two connective regimes, 
to Inset: magnification highlighting the contact peak. 




FIG. 7: Radial distribution function between clusters centers of mass for T = 0.10 and various studied <j> m the cluster phase 
region. Symbols are simulations results, lines are theoretical calculations based on HNC solution of the OZ equation for a 
renormalized Yukawa interactions. Insets: (a) renormalized amplitude A e ff of the Yukawa cluster-cluster interactions and (b) 
clusters number density (in units of a" -3 ) dependence on (j). Dashed lines are linear fits. 



Hence, we can describe the system as composed essentially of repulsively-interacting clusters up to = 0.08, with a 
renormalized amplitude which results from the average of all amplitudes which characterize clusters of different sizes 
and shapes. This result is central for the present paper, and confirms the conjecture put forward in |24j . which was 
then combined with the use of MCT to provide evidence of the existence of a Wigner glass at low cp. 

The picture of interacting Yukawa clusters breaks for cf> > 0.08. The ^^'(r) can no longer be fitted with a Yukawa 
analogue. A careful look at the configurations shown in Fig. [2] shows that, at this point, clusters have completely 
lost their spherical-shape. The radial distribution function shows for tfi > 0.1 a reversal of trend in the main peak 
position and amplitude, which then persists at larger (f>. Moreover, a small peak develops with increasing <j) for 
1 < r/er < 2, indicating the occurrence of some cluster branching events. The Yukawa effective interaction naturally 
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ceases to work in this regime, and a competition between repulsion and increased packing, favoring sometimes the 
formation of intercluster bonds, emerges. This mechanism, which can be considered essentially absent for <j) < 0.08, 
becomes important in this intermediate regime 0.08 < <fi < 0.14, and finally dominant above percolation. These results 
highlight the mechanism by which a crossover between a cluster phase and a percolating state is realized, based on a 
change between repulsion-dominated clusters to a branching regime of clusters due to increased packing. 



V. DYNAMICS: MSD AND ISO-DIFFUSIVITY LINES 



Next, we monitor the dynamics of the system, to find out whether and how we approach dynamic arrest at low T, 
both in the cluster region and in the percolating one. We start by calculating the particle mean squared displacement 
(MSD) and reporting its behaviour in Fig. [8] for T = 0.1 and all studied <fi. We observe a gradual decrease of the 
MSD with increasing (f>. From the time dependence of the MSD one can identify distinct regions. For very short 
times, one observes ballistic motion, followed by a slowing down at short times, which takes place approximately in 
correspondence with the length scale of the attractive bond distance, defined as r m .; n /cr — (2 1 / 100 — 1) ~ 0.007, i.e. 
when two particles are in the minimum of the potential well. This attractive localization by neighbouring particles 
is active for roughly a decade, even at small 4>, in the cluster phase region. It is quite remarkable to find particles 
rattling inside the narrow bonds for such a long time. For intermediate times, particles are able to escape (on average) 
from the bond, as signalled by the fact that they overcome the distance corresponding to the maximum bond distance 
r max /a = (rt — 1) ~ 0.072. Hence, particles explore their neighbourhoods, and a second slowing down emerges. 
This occurs on length scales that can exceed a particle diameter, depending on <f>. These results suggest that the 
motion of a single particle results from the sum of the brownian motion of the whole cluster and of the intra-clustcr 
motion. We notice that only the smallest values of <j> are able to recover a pure diffusive regime at this temperature. 
For larger <f> a sub-diffusive growth of the MSD is found in the time window which can be numerically studied. The 
apparent exponent regulating the subdiffusivity is found to decrease with 0, as expected for an MSD approaching a 
flat plateau. The onset of subdiffusive behaviour signals that arrest is close-by, and located at slightly lower T, both 
for the cluster phase region and for the percolating one. In both cases, particles are trapped at first by the attractive 
bonds with neighbouring particles and secondly by a larger localization length. In the case of a cluster phase, such 
length provides the distance where clusters can rattle (see Movies provided as supplementary information [38 ). For 
the percolating states, the localization length does not decrease much and the arrested state can be identified as a 
gel[4]. Interestingly, the transition between the two arrested states is continuous from the MSD point of view. 




FIG. 8: MSD < r 2 > vs time for T = 0.1 and various <j>. From top to bottom, (f> = 
0.005,0.02,0.04,0.08,0.10,0.125,0.14,0.16,0.20. The dashed line indicates diffusive behaviour. Horizontal lines refer to the 
squared minimum r^ in and maximum bond distance r^ax (see text). 

To further prove that arrest takes place also in the non-percolating regime, we show in Fig.[9]the MSD for = 0.125 
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and various studied T. We remark that this is an isochore along which the system never percolates, at all studied T, 
but is the closest, among the studied <fi, to the percolation transition. Similar results are found at lower <p. At very 
low T, the MSD approaches a flat behaviour at long times, so that true arrest occurs. The low-T curves are calculated 
within a time- window where the energy slowly drifts, and very slow aging effects are thus present. These effects can 
only slow further the relaxation, so that the reported curves can be considered an overestimation of the true MSDs. 
The localization length for the arrested state at this tfi can be estimated around « 0.8<r. The effect of temperature is 
also visible in the amplitude of the vibration within the bond, which decreases with T, as seen in the height of the 
inflection around r m i n . 

We note in passing that, in principle, we can also calculate the MSD of clusters (not shown), by monitoring their 
center of mass trajectory. Due to the fact that clusters undergo breaking and reforming events during the simulation 
time, statistics is poor. However, in all studied cases, the average MSD of clusters center of mass is always found 
below that of a single particle. This provides a hint of the fact that, in the explored time window, most of the 
derealization is provided by rotational motion of the isolated clusters (sec also Movies provided in the supplementary 
information |38j ) . 




FIG. 9: MSD < r 2 > vs time for 4> = 0.125 and various T. From top to bottom, T = 0.2, 0.15, 0.125, 0.10, 0.07, 0.05. Horizontal 
lines refer to the squared minimum r^ in and maximum bond distance r^ aa , (see text). 

From the long-time behaviour of the MSD, we can extract the self-diffusion coefficient, defined as D = 
lim t ^ 00 (r 2 )/6t, for those state points where a diffusive long-time regime can be clearly identified. We can draw 
iso-diffusivity lines, i.e. loci in the phase diagram with constant diffusion coefficient, expressed as a fraction of the 



bare diffusion coefficient Dq. We plot some iso-diffusivity lines in Fig. 10 together with the phase diagram reported 
above. We clearly find that the iso-D lines follow at small <f> the shape of the cluster phase boundary, while at larger <j> 
they follow the percolation boundary. From previous studies [351 HOI SU H2] , we know that the shape of iso-diffusivity 
lines does not change much with approaching distance to the D = 0-linc, which can be identified with the ideal arrest 
transition. Hence, if one could extrapolate, the arrest line would be somehow parallel to the iso-D lines, signaling 
that arrest is mainly temperature-driven. Again, we notice the continuous shape of the iso-diffusivity lines across the 
percolation transition. This allows us to identify low-T states as arrested, or approaching arrest at low T, indepen- 
dently of the presence of a percolating network, the latter condition being discriminant for determining the nature of 
the two arrested states. Hence, we can conclude that also the non-percolating states undergo dynamical arrest at low 
T in full analogy with the percolating ones. Thanks to the identification of the dominant cluster-cluster interactions 
discussed above, arrest at low <j> can therefore be interpreted as a Wigner glass of clusters, and it is clearly distinct 
from arrest driven by the formation of a spanning cluster (gel) . 

Due to the presence of subdiffusion and of aging effects at low T, we can only probe a limited window in D/Do. The 
behaviour of D vs T (not shown) displays a rapid decrease below T — 0.20, then followed by a slower (non-Arrhenius) 
decrease for T < 0.125. Due to the limited available data and the difficulty to extract D for the low-T states, our 
analysis can not be detailed. However, the apparent non-Arrhenius behaviour might indicate the presence of a finite- 
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FIG. 10: Phase diagram, partially redrawn from Fig. [TJ with added iso-diffusivity lines, for two selected values of D/Dq. 

T arrest line, at least for the larger densities, and in particular in the gel regime. This is different from what was 
reported for other gel-forming systems, based on patchy or limited-valence interactions, where Arrhenius behaviour 
has been observed EHl H5] . compatible with an ideal gel state occuring only at T — > 0(11 HB]. Indeed, in the patchy 
particle models, the relaxation and approach to a gel state are essentially controlled by the single-bond lifetime [42 j 
and bonds were shown to be independent. For the present case, the interplay between different types of interactions, 
and in particular their long-range nature, might be responsible for the fact that collective bond rearrangements are 
needed in order to restructure the system. 

VI. CONCLUSIONS 

In this manuscript, we have reported extensive Brownian Dynamics simulations of a model potential suitable 
to describe the interactions between charged colloidal particles in apolar solvent, in the presence of an additional 
short-range depletion attraction. The long-range nature of the repulsive screened Coulomb interactions is taken into 
account numerically by performing Ewald summation. We have studied a wide region of the (<j),T) plane, covering 
the evolution of the system from a fluid of monomers to a percolating network. In between these two limits, a region 
of stable finite-size clusters is found. These clusters are spherical at low 4>, acquiring more and more elongated shape 
with increasing <j>, due to the role of cluster-cluster interactions. 

We have provided evidence that particles inside clusters can acquire local order, but no long-range order, both inside 
clusters and among different clusters, is present. Most importantly, we have been able to quantify the cluster-cluster 
interactions in terms of an effective Yukawa potential, with the same screening length as for particle-particle interac- 
tions, and a renormalized amplitude, confirming a previously hyphothesized scenario |24j . This happens independently 
of the clusters polydispersity in size and shape. However, when <fi increases, the situation changes, and cluster-cluster 
interactions become more subtle, due to the increasing number of branching events and to the elongation of the 
clusters. 

At low temperatures, we find evidence of dynamic arrest, by monitoring the MSD and the particles self-diffusion 
coefficient. The arrest mechanism appears to be continuous across all studied <f> and is mostly driven by temperature. 
However, the distinct nature of the two different low-T states, namely finite-size clusters or spanning network, allows 
us to unambiguously identify the presence of two different non-ergodic states in these colloidal systems. A Wigner glass 
of clusters exists at low/intermediate cf>, stabilized by the renormalized Yukawa cluster-cluster interactions discussed 
above, while a gel state, stabilized by the presence of a percolating, long-lived network takes place at larger </>. It will 
be interesting in future studies to compute the viscoelastic response of these two different non-ergodic systems, that 
should manifest extremely different rheological properties. The accurate knowledge of the phase diagram, provided in 
this work, will allow us to choose the best conditions where such response is measurable and to compare it to those 
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of other gel or glass-forming systems, as well as to experimental results. 
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